#ifndef __FIRSTORDER_H__
#define __FIRSTORDER_H__
/// Includes projects
#include "reduce.h"
#include "memoryRoutines.h"
#include "svmKernels.h"
/// Includes system
#include <cfloat>

__device__ unsigned int count = 0;
__shared__ bool isLastBlockDone;

#define TAU 1e-12f
#define NONE 0x0
#define LOW 0x1
#define UP 0x2

template<bool iLowCompute, bool iHighCompute, class Kernel>
__global__ void	firstOrder(//float* devData, int devDataPitchInFloats,
						   float* devTransposedData, int devTransposedDataPitchInFloats,
						   float* devLabels,
						   int nPoints, int nDimension,
						   float* devAlpha,
						   float* devF,
						   float sAlpha1Diff, float sAlpha2Diff,
						   int iLow, int iHigh,
						   float parameterA, float parameterB, float parameterC,
						   float* devCache, int devCachePitchInFloats,
						   int iLowCacheIndex, int iHighCacheIndex,
						   int* devLocalIndicesRL, int* devLocalIndicesRH,
						   float* devLocalFsRL, float* devLocalFsRH,
						   float* devKernelDiag,
						   float* devSelfDot,
						   float* devResult,
						   float cost,
						   Kernel kf)
{
	extern __shared__ float xIHigh[];
	float * xILow;
	unsigned int tid = threadIdx.x;
	unsigned int bid = blockIdx.x;
	unsigned int globalIndex = __mul24(blockDim.x, bid) + tid;
	
	float f;
	int reduceFlag = NONE;
	__shared__ int tempLocalIndices[BLOCKSIZE];
	float * xILow;
	unsigned int tid = threadIdx.x;
	unsigned int bid = blockIdx.x;
	unsigned int globalIndex = __mul24(blockDim.x, bid) + tid;
	
	float f;
	int reduceFlag = NONE;

	{
		if (iHighCompute) { xILow = &xIHigh[nDimension]; } 
		else { xILow = xIHigh; }

		if (iHighCompute) {
			scatterRowVector(tid, devTransposedData+ devTransposedDataPitchInFloats*iHigh + tid, nDimension, xIHigh);
			
			//scatterTextureRowVector(tid, xIHigh, iHigh, nDimension);
			//__syncthreads();
		}
		if (iLowCompute) {
			scatterRowVector(tid, devTransposedData+ devTransposedDataPitchInFloats*iLow + tid, nDimension, xILow);

			//scatterTextureRowVector(tid, xILow, iLow, nDimension);
			//__syncthreads();
		}
		if(iHighCompute || iLowCompute) {
			__syncthreads();
		}
	}

	if (globalIndex < nPoints) {
		{
			float alpha = devAlpha[globalIndex];
			float label = devLabels[globalIndex];

			if(alpha > 0.0f && alpha < cost) {
				reduceFlag = LOW | UP;//I0:SV
			} else if(alpha == 0.0f) {
				if( label < 0.0f ) reduceFlag = LOW;//I4 negative nonSV
				else reduceFlag = UP;//I1 positive nonSV
			} else { //alpha == cost
				if( label < 0.0f) reduceFlag = UP;//I2 negative bounded SV
				else reduceFlag = LOW;//I3 positive bounded SV
			}
		}

		{
			//if (!iHighCompute) {highKernel = devCache[(devCachePitchInFloats * iHighCacheIndex) + globalIndex];}
			//if (!iLowCompute) {lowKernel = devCache[(devCachePitchInFloats * iLowCacheIndex) + globalIndex];}

			float highKernel;
			float lowKernel;

			//if(dynamic_cast<Gaussian *>(&kf)) {
			//	if(globalIndex == 0)
			//		printf("match\n");
			//}

			//if(typeid(&kf) != typeid(Gaussian*)) {
			//	//if(globalIndex == 0)
			//	//	printf("match\n");
			//}

			f = devF[globalIndex];
			if (iHighCompute && iLowCompute) {
				kf.dualK(globalIndex, xIHigh, xILow, parameterA, parameterB, parameterC, highKernel, lowKernel, devSelfDot, iHigh, iLow);
				//kf.dualKernel(devData + globalIndex, devDataPitchInFloats, devData + (devDataPitchInFloats * nDimension), xIHigh, 1, xILow, 1, parameterA, parameterB, parameterC, highKernel, lowKernel);
			} else if (iHighCompute && !iLowCompute) {
				lowKernel = devCache[(devCachePitchInFloats * iLowCacheIndex) + globalIndex];
				highKernel = kf.K(globalIndex, xIHigh, parameterA, parameterB, parameterC, devSelfDot, iHigh);
				//highKernel = kf.kernel(devData + globalIndex, devDataPitchInFloats, devData + (devDataPitchInFloats * nDimension), xIHigh, 1, parameterA, parameterB, parameterC);
			} else if (iLowCompute && !iHighCompute) {
				highKernel = devCache[(devCachePitchInFloats * iHighCacheIndex) + globalIndex];
				lowKernel = kf.K(globalIndex, xILow, parameterA, parameterB, parameterC, devSelfDot, iLow);
				//lowKernel = kf.kernel(devData + globalIndex, devDataPitchInFloats, devData + (devDataPitchInFloats * nDimension), xILow, 1, parameterA, parameterB, parameterC);
			} else {
				highKernel = devCache[(devCachePitchInFloats * iHighCacheIndex) + globalIndex];
				lowKernel = devCache[(devCachePitchInFloats * iLowCacheIndex) + globalIndex];
			}

			f += sAlpha1Diff * highKernel;
			f += sAlpha2Diff * lowKernel;
			
			if (iLowCompute) devCache[(devCachePitchInFloats * iLowCacheIndex) + globalIndex] = lowKernel;
			if (iHighCompute) devCache[(devCachePitchInFloats * iHighCacheIndex) + globalIndex] = highKernel;
		}
		devF[globalIndex] = f;
	}
	
	if ((reduceFlag & LOW)==NONE) {
		tempLocalFs[tid] = NINFTY; //Ignore me
	} else {
		tempLocalFs[tid] = f;
		tempLocalIndices[tid] = globalIndex;
	}
	__syncthreads();

	argmaxReduce<BLOCKSIZE>(tempLocalFs, tempLocalIndices, tid);
	if (tid == 0) {
		devLocalIndicesRL[bid] = tempLocalIndices[0];
		devLocalFsRL[bid] = tempLocalFs[0];
	}
	__syncthreads();

	if ((reduceFlag & UP)==NONE) {
		tempLocalFs[tid] = INFTY; //Ignore me
	} else {
		tempLocalFs[tid] = f;
		tempLocalIndices[tid] = globalIndex;
	}
	__syncthreads();

	argminReduce<BLOCKSIZE>(tempLocalFs, tempLocalIndices, tid);
	if (tid == 0) {
		devLocalIndicesRH[bid] = tempLocalIndices[0];
		devLocalFsRH[bid] = tempLocalFs[0];
		
		__threadfence();
		int value = atomicAdd((int*)&count,1);	//declare "value" prevent race condition
		isLastBlockDone = (value == (gridDim.x-1));
	}
	__syncthreads();

	if(isLastBlockDone) {
		__shared__ int iHigh,iLow;
		__shared__ float bHigh,bLow;

		if (tid < gridDim.x) {
			tempLocalIndices[tid] = devLocalIndicesRH[tid];
			tempLocalFs[tid] = devLocalFsRH[tid];
		} else {
			tempLocalFs[tid] = INFTY;
		}
		for (int i = tid + blockDim.x; i < gridDim.x; i += blockDim.x) {
			if(devLocalFsRH[i] < tempLocalFs[tid]) {
				tempLocalFs[tid] = devLocalFsRH[i];
				tempLocalIndices[tid] = devLocalIndicesRH[i];
			}
		}
		__syncthreads();
		argminReduce<BLOCKSIZE>(tempLocalFs, tempLocalIndices, tid);
		if(tid == 0) {
			iHigh = tempLocalIndices[0];
			bHigh = tempLocalFs[0];
		}
		__syncthreads();

		if (tid < gridDim.x) {
			tempLocalIndices[tid] = devLocalIndicesRL[tid];
			tempLocalFs[tid] = devLocalFsRL[tid];
		} else {
			tempLocalFs[tid] = NINFTY;
		}
		for (int i = tid + blockDim.x; i < gridDim.x; i += blockDim.x) {
			if(tempLocalFs[tid] < devLocalFsRL[i]) {
				tempLocalIndices[tid] = devLocalIndicesRL[i];
				tempLocalFs[tid] = devLocalFsRL[i];
			}
		}
		__syncthreads();
		argmaxReduce<BLOCKSIZE>(tempLocalFs, tempLocalIndices, tid);
		if(tid == 0) {
			iLow = tempLocalIndices[0];
			bLow = tempLocalFs[0];
		}
		__syncthreads();

		{
			float* highPointer = devTransposedData + (iHigh * devTransposedDataPitchInFloats);
			float* lowPointer = devTransposedData + (iLow * devTransposedDataPitchInFloats);
			kf.parallelKernel(tid, highPointer, highPointer + nDimension, lowPointer, tempLocalFs, parameterA, parameterB, parameterC);
		}

		if (tid == 0) {
			count = 0;
			float eta = devKernelDiag[iHigh] + devKernelDiag[iLow];
			float kernelEval = tempLocalFs[0];

			eta = eta - 2.0f*kernelEval;
			if (eta <= 0.0f) eta = TAU;

			float alpha1Old = devAlpha[iHigh];
			float alpha2Old = devAlpha[iLow];
			float lowLabel = devLabels[iLow];
			float highLabel = devLabels[iHigh];
			float delta = lowLabel*(bHigh-bLow)/eta;
			//float sign = devLabels[iHigh] * lowLabel; //if sign < 0 => y1 != y2
			float alpha1New;
			float alpha2New;
			
			if (lowLabel * highLabel < 0.0f) {
				float alphaDiff = alpha2Old - alpha1Old;
				alpha1New = alpha1Old + delta;
				alpha2New = alpha2Old + delta;
				if(alphaDiff < 0.0f) {
					if(alpha2New < 0.0f) {
						alpha2New = 0.0f;
						alpha1New = -alphaDiff;
					}
					if(alpha1New > cost) {
						alpha1New = cost;
						alpha2New = cost + alphaDiff;
					}
				} else {
					if(alpha1New < 0.0f) {
						alpha1New = 0.0f;
						alpha2New = alphaDiff;
					}
					if(alpha2New > cost) {
						alpha2New = cost;
						alpha1New = cost - alphaDiff;
					}
				}
			} else {
				alpha1New = alpha1Old - delta;
				alpha2New = alpha2Old + delta;
				float alphaSum = alpha2Old + alpha1Old;
				if(alphaSum > cost) {
					if(alpha1New > cost) {
						alpha1New = cost;
						alpha2New = alphaSum - cost;
					}
					if(alpha2New > cost) {
						alpha2New = cost;
						alpha1New = alphaSum - cost;
					}
				} else {
					if(alpha2New < 0.0f) {
						alpha2New = 0.0f;
						alpha1New = alphaSum;
					}
					if(alpha1New < 0.0f) {
						alpha1New = 0.0f;
						alpha2New = alphaSum;
					}
				}
			}
			devAlpha[iLow] = alpha2New;
			devAlpha[iHigh] = alpha1New;

			//devResult[0] = alpha2Old;
			//devResult[1] = alpha1Old;
			devResult[0] = lowLabel*(alpha2New - alpha2Old);
			devResult[1] = highLabel*(alpha1New - alpha1Old);
			devResult[2] = bLow;
			devResult[3] = bHigh;
			//devResult[4] = alpha2New;
			//devResult[5] = alpha1New;
			devResult[4] = __int_as_float(iLow);
			devResult[5] = __int_as_float(iHigh);
		}
	}
}


//#include "firstOrder.cu"

#endif
